//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   Quadrature.cc
 * \author Jeremy Roberts
 * \date   06/19/2011
 * \brief  Quadrature member definitions.
 * \note   Copyright (C) 2011 Jeremy Roberts. Derived from Denovo.
 */
//---------------------------------------------------------------------------//
// $Rev:: 141                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-09-14 19:28:23 +0000 (Wed, 14 Sep 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <cstdio>

#include "angle/Quadrature.hh"

namespace slabtran
{

//---------------------------------------------------------------------------//
// CONSTRUCTOR
//---------------------------------------------------------------------------//
/*!
 * \brief Constructor.
 */
template <class D>
Quadrature<D>::Quadrature(int    sn_order,
                          double norm)
    : d_sn_order(sn_order)
    , d_norm(norm)
{
    Require (d_sn_order % 2 == 0);
    Require (d_norm > 0.0);
}

// Angle Accessor general definitions

//! Return the angles.
template <class D>
const Vec_Dbl& Quadrature<D>::angle(int angle) const
{
    Require (angle >= 0 && angle < d_angles.size());
    return d_angles[angle];
}

//! Return the angle for a given ordinate.
template <class D>
double Quadrature<D>::angle(int angle, int m) const
{
    Require (angle >= 0 && angle < d_angles.size());
    Require (m >= 0 && m < d_angles[angle].size());
    return d_angles[angle][m];
}

//! Return the dimensionality.
template <class D>
int Quadrature<D>::dimensionality() const { return d_angles.size(); }

// Angle accessor specializations for 1D

//! Return the angles (a reference to the vector)
template <>
const Vec_Dbl& Quadrature<_1D>::angle(int angle) const
{
    Require (angle==MU);
    return d_angles;
}

//! Return the angle for a given ordinate.
template <>
double Quadrature<_1D>::angle(int angle, int m) const
{
    Require (angle==0);
    Require (m >= 0 && m < d_angles.size());
    return d_angles[m];
}

//! Return the dimensionality.
template<>
int Quadrature<_1D>::dimensionality() const { return 1; }


//---------------------------------------------------------------------------//
/*!
 * \brief Virtual destructor.
 */
template <class D>
Quadrature<D>::~Quadrature()
{
}

//---------------------------------------------------------------------------//
/*!
 * \brief Display the quadrature.
 *
 * \todo The quadratures should print in the same way they are stored and
 *       and accessed.
 */
template <class D>
void Quadrature<D>::display() const
{
	using std::cout;
	using std::endl;
	using std::printf;
	Require ( d_angles.size() == _3D::dimension );
	cout << endl;
	name(); // print name of quad set
	cout << endl << "The quadrature directions and weights are:"
		 << endl << endl;
	cout << "   m         mu              eta             xi              wt       " << endl;
	cout << "  ---   --------------  --------------  --------------  --------------" << endl;
	double sum_wt = 0.0;
	for ( int ix = 0; ix < d_angles[MU].size(); ++ix )
	{
		printf ("%4i    %13.10f   %13.10f   %13.10f   %13.10f   \n",
				ix, d_angles[MU][ix], d_angles[ETA][ix], d_angles[XI][ix], d_wt[ix] );
		sum_wt += d_wt[ix];
	}
	cout << endl << "  The sum of the weights is " << sum_wt << endl;
	cout << endl;
}
template <>
void Quadrature<_2D>::display() const
{
	using std::cout;
	using std::endl;
	using std::printf;
	Require ( d_angles.size() == _2D::dimension );
	name(); // print name of quad set
	cout << endl << "The quadrature directions and weights are:"
			<< endl << endl;

	cout << "   m         mu              eta             wt       " << endl;
	cout << "  ---   --------------  --------------  --------------" << endl;

	double sum_wt = 0.0;
	for ( int ix = 0; ix < d_angles[MU].size(); ++ix )
	{

		printf ("%4i    %13.10f   %13.10f   %13.10f   \n",
				ix, d_angles[MU][ix], d_angles[ETA][ix], d_wt[ix] );

		sum_wt += d_wt[ix];
	}
	cout << endl << "  The sum of the weights is " << sum_wt << endl;
	cout << endl;
}
template <>
void Quadrature<_1D>::display() const
{
	using std::cout;
	using std::endl;
	using std::printf;
	Require ( d_angles.size() == _1D::dimension );
	name(); // print name of quad set
	cout << endl << "The quadrature directions and weights are:"
			<< endl << endl;
	cout << "   m         mu              wt       " << endl;
	cout << "  ---   --------------  --------------" << endl;
	double sum_wt = 0.0;
	for ( int i = 0; i < d_angles.size(); ++i )
	{
		printf ("%4i    %13.10f   %13.10f \n", i, d_angles[i], d_wt[i]);
		sum_wt += d_wt[i];
	}
	cout << endl << "  The sum of the weights is " << sum_wt << endl;
	cout << endl;
}
//! instantiations
template class Quadrature<_1D>;
template class Quadrature<_2D>;
template class Quadrature<_3D>;



} // end namespace slabtran

//---------------------------------------------------------------------------//
//                 end of Quadrature.cc
//---------------------------------------------------------------------------//
